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We present a method to perform stability analysis of nonequilibrium fixed points appearing in 
self-consistent electron transport calculations. The nonequilibrium fixed points are given by the 
self-consistent solution of stationary, nonlinear kinetic equation for single-particle density matrix. 
We obtain the stability matrix by linearizing the kinetic equation around the fixed points and 
analyze the real part of its spectrum to assess the asymptotic time behavior of the fixed points. 
We derive expressions for the stability matrices within Hartree-Fock and linear response adiabatic 
time-dependent density functional theory. The stability analysis of multiple fixed points is performed 
within the nonequilibrium Hartree-Fock approximation for the electron transport through a molecule 
with a spin-degenerate single level with local Coulomb interaction. 
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I. INTRODUCTION 

The existence of nonunique steady state for nonequilib- 
rium systems of correlated quantum or classical particles 
is an interesting and open fundamental problem. 0, 0] 
Multiple steady states in nano junctions lead to bista- 
bilities and hysteresis loops in the current- voltage char- 
acteristics and this currently is a much debated the- 
oretical issue. [3-5] Density functional theory (DFT) 
and Hartree-Fock based nonequilibrium Green's func- 
tions (NEGF) electron transport calculations, which 
are widely used nowadays, solve nonlinear system of 
equations for nonequilibrium electron density via self- 
consistent iterations. Such kind of nonlinear prob- 
lems may have multiple solutions. [H, [T^ In equilibrium 
case the correct physical solution corresponds to a min- 
imum of the ground state energy. The situation is less 
clear in nonequilibrium where the system is open and the 
minimum energy arguments are not applicable anymore. 
In this paper, we discuss the appearance of multiple fixed 
points and present a method to eliminate unphysical solu- 
tions of steady-state nonequilibrium self-consistent elec- 
tron transport problem. 

Let us consider a finite quantum system (e.g., a 
molecule) connected to two macroscopic particle reser- 
voirs or thermal baths (e.g., metal electrodes). By pro- 
jecting out bath degrees of freedom we obtain the kinetic 
equation for the reduced density matrix p(t) of the em- 
bedded system 



d , , 



Cp(t), 



(i) 



where C is a non-Hermitian Liouvillian. We assume 
that the dynamics is markovian and the system is au- 
tonomous, i.e. C does not depend explicitly on time. In 
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this paper we will focus on nonequilibrium steady state. 
Like an equilibrium represents stationary state of a closed 
system, a nonequilibrium steady state is the stable, time- 
invariant state of an open system. If the dynamics gen- 
erated by C is linear, then a nonequilibrium steady state 
can be unambiguously defined as a state when the left 
side of the equation (QJ becomes zero. However, often 
times for practical calculations we involve the mean-field 
approximation (Hartree-Fock or DFT) and the Liouvil- 
lian becomes a functional of the reduced density matrix 



C = C[p], 



(2) 



Therefore the kinetic equation ([T]) becomes nonlinear and 
the issue of stability of the solution becomes pivotal. 

Let us give a few formal definitions from the theory of 
dynamical dissipative systems. (Til Il2j which are relevant 
to electron transport problem. The density matrices ~p at 
which 



C[p]p = 0, 



(3) 



are called the fixed points of the system. Since Eq.Q 
is nonlinear, it generally has multiple solutions, which 
may or may not be steady state (i.e. stable fixed point). 
To understand whether or not the fixed point is stable 
wc need to perform stability analysis commonly used for 
dynamical nonlinear systems. We expand the density 
around the fixed point 



p(t) = p + Sp(t) 
and linearize the Liouville equation 



where 



-Sp(t) =A6p(t), 



dp 



p=p 



(4) 



(5) 



(6) 



is the stability matrix. If A is a so-called Hurwitz ma- 
trix, i.e., if all eigenvalues Aj of A satisfies the conditions 
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Re(Ai) < 0, then the fixed point ~p is asymptotically sta- 
ble and it is the true steady state of the system. If at 
least one eigenvalue has positive real part, the solution 
is unstable and can not be a steady state, since even an 
infinitesimally small variation of the fixed point density 
matrix drives the solution away from the fixed point ex- 
ponentially in time. In our paper we demonstrate that 
there are several fixed points and multiple steady states 
in some rather typical cases of electronic transport cal- 
culations. 

The rest of the paper is organized as follows. In Section 
II, we discuss nonequilibrium fixed points, linearization 
and stability matrix for self-consistent electron transport 
problem and obtain the general expression for the stabil- 
ity matrix. In section III, we apply the method to out 
of equilibrium Anderson model and demonstrates that it 
has multiple (stable and unstable) fixed points. Conclu- 
sions are given in Section IV. In Appendix A, we derive 
the kinetic equation for the reduced density matrix of the 
embedded system. The adiabatic time-dependent Kohn- 
Sham DFT expressions for stationary Fock matrix and 
its fluctuating part are given in ppendix B. We use natu- 
ral units throughout the paper: h = Ub = |e| = 1, where 
— |e| is the electron charge. 



II. NONEQUILIBRIUM FIXED POINTS, 
LINEARIZATION AND STABILITY MATRIX 

Let us consider a molecule connected to two electrodes. 
We partition the system into five parts: the molecule it- 
self, the left /right macroscopically large leads (environ- 
ment), and the left/right finite buffer zones between the 
molecule and the environment. The Hamiltonian is writ- 
ten in the following form: 

T~L = Hm + He + Hb + Hmb + Heb ■ (7) 

The environment and the buffer zones are described by 
the noninteracting Hamiltonians 

H E = E £kal a a k a, (8) 

a,k£L,R 



H B = E eba\ a a b< y. (9) 

a.beL.R 

Here e k denote the continuum single-particle spectra of 
the left (fc 6 L) and right (fc € R) lead states, a\ a (a*^) 
create (annihilate) electron with spin a =~f, ], in the lead 
state k. The buffer zones have discrete energy spectrum 
Eb with corresponding creation and annihilation opera- 
tors al a and a ba . 

The molecular Hamiltonian is taken in the most gen- 
eral form and contains electronic kinetic energy, electron- 
ion interaction, and Coulomb interaction between elec- 



trons: 

H M = ^2,Tija\ a a ja + 7; E ^^\ mn ) a \a a Lc,' a na'a jrj , 

ija ijmn a a' 

(10) 

where 

T nm = J rfr^„(r)(-iv 2 + V c _i(r))0 m (r) (11) 

and 

(ij\mn) = J drdr>i(r)^(r) j ^ m (r>»(O (12) 

are matrix elements computed in some orthogonal real 
basis (4>i\4>j) — Sij. Here a\ a and dj ff are creation and 
annihilation operators for electron with spin a in molec- 
ular state \4>i)- Hereinafter, the index b refers to dis- 
crete single-particle states in either left or right buffer 
zones, whereas the indices i,j,m,n represent states in 
the molecular space. 

The buffer-environment and molecule-buffer coupling 
have the standard tunneling form: 

Heb = E ( v bkal a a k(7 +h.c.)+ ^ {vbka\ a a ka+h.c), 

a,bk£L a,bkeR 

(13) 

H M b = E (tibaOfoOur + h.c.). (14) 

icr,b£L,R 

The Liouville equation for the total density matrix x{t) 

is: 

ix(t) = [H,x(t)}. (15) 

If we project out the environment degrees of freedom and 
make the standard assumptions (Born-Markov and rotat- 
ing wave approximations [13J - the details of the deriva- 
tion are shown in Appendix A) we get the following ki- 
netic equation 

ip(t) = [H,p(t)] + Up(t) (16) 

for the reduced density matrix for the embedded molecule 
(i.e. for the molecule and the buffer zones) 

p(t)=Ti EX (t). (17) 

Here the Hamiltonian H includes the Lamb shift of the 
single-particle levels of the buffer zones 

H = H M + H M b + E( £b + &b)al a a ba , (18) 

a,b 

and the non-Hermitian dissipator is given by standard 
Lindblad form 

= E E ( 2L b^pml, - {L\^L b ^, P {t)}) 

a,b /i=l,2 

(19) 
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with the following Lindblad operators 

Lbai = Vlb(l - f b )a ba , L ha2 = V7b/f> a L- (20) 

Here Af, and "f b are real and imaginary parts of the stan- 
dard environment self energy ^ fc \v b k\ 2 /(s b — £fc + i0 + ) 
and f beL/R = [1 + e^/fi^-^/fi)]-!. 

We emphasize that our kinetic equation (IT6)) does 
not employ the second order perturbative treatment of 
molecule-electrode coupling (fT4| . but rather it is the 
second-order in terms of the coupling between the buffer 
zone and the environment (|13[) . We have recently demon- 
strated that for the steady state electron transport cal- 
culations the kinetic equation (ITB1) can be made as ac- 
curate and exact as necessary in practical calculations 
by increasing the density of sin gle- particle buffer states 
b included into the Hamiltonian \l4, [H[ The similar idea 
of the buffer zone between the molecule and the environ- 
ment has been recently proposed in dynamical simula- 
tions of inelastic electron transport. fl6| 

Let us now consider time-evolution of the expectation 
value of an arbitrary operator O: 

(0} t =Tr[p(t)0}. (21) 

Using Eq. (fT5|) we obtain 

j t {0) t = -i{[0,H]) t 
+ E E ( 2 ( L l, 0L ^)t - ({Ll,L b ^O}) t ). (22) 

ab n=l,2 

I 



with Eb = Sb — ijb- Here we include the Lamb shift At, 
into single-particle energy Sb- These equations of motion 
are nonlinear because the Fock matrix Ffj (t) depends on 
the density matrix (t) . 

Setting the left side of equations (f27|) to zero, we obtain 
the stationary Hartree-Fock equations for nonequilibrium 
fixed points. These fixed points correspond to stationary 
single-particle densities, Py, P bi ,P ib , and P bb i, which 
may or may not be steady state densities. To deter- 



We apply the time-dependent Hartree-Fock approxima- 
tion to the two-particle interaction in the molecular 
Hamiltonian 

#m=E^KW- (23) 

crij 

Here F° m (t) is the time-dependent Fock matrix 

KnM = T„ m + E [("Hv W*) - (jii\mj)PtM , 

which depends on the single-particle density matrices 
P°{t) = (a] a a. la ) t , P t3 {t) = J2 P^(f)- (25) 

IT 

The corresponding expression for the Fock matrix in adi- 
abatic time-dependent Kohn-Sham DFT is given in Ap- 
pendix B. The Lindblad equation (|16[) becomes the time- 
dependent Hartree-Fock equation for the nonequilibrium 
open quantum system, when the full many-body molec- 
ular Hamiltonian is approximated by (|2"3")l . We in- 
troduce two additional single-particle density matrices: 

Phi® = <SW>* > Pw(t) = (°t**b<r)t (26) 

and P? b {t) = (Pjg (*))*. By means of (J22J we get the close 
set of time evolution equations for single-particle density 
matrices 



(27) 



I 

mine if these fixed points are asymptotically stable, we 
linearize Eq. (1271) around each fixed point. Substituting 
(here indices a, j3 run over molecular and buffer single 
particle states) 

W)=Kfi+SP^)(t) (28) 

into Eq. (|27[) and retaining only the terms linear in 
6PZ p (t) we get 



n b 

ij/m = ^(o-E^w^-w-E^^w+E^^W' 

3 b ' j 

iitw® = ( E b ~ K)Pw® - E[^'- P ^W - tl ba Pg,{t)]+2i5 w f blb 
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Y^i^JP^t) - F a nj 5P? n {t) + 5F° n {t)P a nj - SF^(t)P- n ] +J2[tib*6Pb j (t) - fjbJPm)] , 

n b 

E b 5Pm -^[F^SP^t) + 5F°(t)P a bj ] -Y,**a>JPvb(t) + ^t* jba SI^(t), 

3 V 3 

(E b - E* b ,)SP b ° b , (t) - [UvaSPm - tlJPg, (t)] . (29) 

i 

i 



Here the Fock matrix at the fixed point is 



(30) 



and its time-dependent fluctuation around this fixed 
point is 

SFZ m (t) = ]T [(nm|ij)M> tf (t) - (m\mj)5P°(t)] . (31) 



The adiabatic time-dependent Kohn-Sham DFT expres- 
sions for stationary Fock matrix and its fluctuating part 
are given in appendix B. 

The system of equations (f2U)l can be rewritten as the 
set of linear differential equation 



d 
dt 



SP*p(t) = E E A $, a>r&P*p>{t)- (32) 

at' /3' a' 



Here A^g a ,g, is the spin-dependent stability matrix. 
Each element of this stability matrix can be readily ob- 
tained from (f2U)l . In general case it is sparse, complex 
and non-Hermitian matrix. Now we need to find eigen- 
values Xi of the stability matrix a ,^, and analyze 
their real parts. If Re(Ai) < for all eigenvalues in the 
spectrum, then the fixed point P a p is asymptotically sta- 
ble as t — > oo and, therefore, it is the true steady state 
of the system. If at least one eigenvalue has a positive 
real part, the solution becomes unstable along this mode 
and it is not a steady state. [Ill G3 The purely imaginary 
eigenvalues of stability matrix correspond to the peri- 
odically oscillating fixed points which may be relevant to 
dynamical picture of Coulomb blockade regime. [TtJ If one 
of Xi becomes zero, the system undergoes zero-eigenvalue 
bifurcation which can be saddle-node, transcritical, or 
pitchfork type bifurcation. fl2j 



III. EXAMPLE CALCULATIONS 

To illustrate the theory we consider electron trans- 
port through a molecule with a spin-degenerate single 
level with local Coulomb interaction (so called Anderson 
model). The Hamiltonian is given by 
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FIG. 1: Schematic illustration of the model system used in 
the electron transport calculations. In the upper part of 
the figure only the first left 6i (1 £ L) and the last right 
eat (N £ R) atoms from the buffer zones are attached to 
the environment. After the diagonalization of the buffer zone 
Hamiltonian, each energy level Sb is connected to the dissipa- 
tors and the molecule (lower part of the figure). 



where the molecular Hamiltonian is 



(34) 



In our calculations left and right buffer zones are 
modeled as a finite chain of N atoms, characterized by 
the hopping parameter Vh and the on-site energy €l,r 
(Fig. [TJ. Thus, the energy spectrum of the each buffer is 
given by 



£b = £l,r + 2V h cos 



■rrb 



N + 1 



b = l,...,N. (35) 



If V c is the spin independent coupling between the 
molecule and the edge buffer site, then the tunneling ma- 
trix elements in Eq. (|33[) are 



H = H M - E ifc ( a ^ ab<T + ^" c ) + E £ f> a L a ^' ( 33 ) t b = V c 



ab 



ah 



2 ( Tib 

■ sin 



N + 1 



N + l 



, b=l,...,N. (36) 



1 
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FIG. 2: Upper panel: The graphical solution of stationary 
nonequilibrium Hartree-Fock equations (I38[) for different val- 
ued of U. Lower panel: the fixed point electron densities 
rife = P sa (k = 1,2,3) as functions of U. Black circles indi- 
cates saddle-node bifurcation points. 



The parameter "/ b in Lindblad master equation is taken to 
be equal to the distance between neighbor energy levels 
in the buffer zones, i.e., ^ b — e b — £b+i- 

In our calculations we use the same parameters as 
in 0. Namely, Vh = —0.5, V c = —0.35, and the molecu- 
lar orbital energy is £0 =0. The buffer on-site energies 
are shifted by the applied voltage bias (V = Vl — Vr): 
£lm = 0.3 + V l ,r, where V L = 1.5 and Vr = 0.0. 
The leads are half-filled so that the Fermi levels of lead 
L,R are positioned at £l.r- The inverse temperature is 
(3 = 90. We use N = 400 and this choice will be justified 
below. Here we keep U as a variable quantity and illus- 
trate how the number and properties of nonequilibrium 
fixed points depend on the strength of electron-electron 
interaction. 

The time-dependent Hartree-Fock dynamics for the 
model Hamiltonian (|33p is fully characterized by the fol- 



lowing (spin independent) single-particle density matrix: 
P sa (t) = (aj T a cr ) t , 

P.b(t) = (aL°*>t> p bb'(t) = (4 ff <«to)i- (37) 

To find the fixed point densities, P a f}, we solve the 
nonequilibrium stationary Hartree-Fock equations (i.e. 
the system of equations (|27p with the time-derivatives 
of the density matrix set to zero): 

J2h (Psb-Pbs) =0, 

b 

t b P ss + (e + UP SS - E* b )P sb - tb'Pvb = 0, 

b' 

t b P sb , - t v P b s - {E b - E* v )P bb > = 2i5 w fblb- (38) 

These Hartree-Fock equations are nonlinear with respect 
to the fixed point electron density in the molecule, n = 

P ss • 

By numerical solution of the Hartree-Fock equations 
we have found that there is a range of the Coulomb in- 
teraction strength parameters, U, for which Eq. (|38[) have 

multipole fixed point solutions P a p (k = 1, 2, . . .). In the 
upper panel of Fig. [5] we show the graphical solution of 
Eq. (f38|) for the different values of U. In this plot the 
crossings of the straight and curved lines give molecular 
densities n = P ss corresponding to different fixed point 
solutions of Eq. (|38[) . We see that three cases are possible 
depending on the value of U: when 1.988 < U < 2.296 
there exist three fixed points; at the ends of the inter- 
val we have two fixed points; and outside the interval 
there is only one fixed point. Hereafter we will number 
fixed point molecular densities in ascending order, i.e., 
ni < n 2 < n 3 . 

In the lower panel of Fig. [5] we show how the molecular 
densities depend on the Coulomb interaction strength. 
The "middle" density ni demonstrates the strong U de- 
pendency and it exists only when 1.988 < U < 2.296. 
The fixed point corresponding to n 2 collides with that 
corresponding to 713 (ni) and annihilate it at the left 
(right) end of this interval for U. This is so-called saddle- 
node bifurcation point [l2| (indicated by black circles in 

Fig.rj. 

Let us understand which of these three fixed point 
solutions are asymptotically stable (i.e., nonequilibrium 
steady states). For this purpose we construct the sta- 
bility matrix (|32p for our model Hamiltonian. The only 
non-zero matrix elements of the stability matrix are 

A S s,sb — ~itbi ^-ss.bs — Ubi ^sb,b'b — Hv 

A sb ,ss = ~i{t b + UP sb ), A sb ,sb = -i[(e + UP SS ) - El], 
Abs.ss = i(h + UP bs ), A bs ,bs = i[(e + UP SS ) - E b ], 

Abs,bb' = —itb'j Abb',sb' = itb,Abb',bs — —itb', 

A WtW = -i(E b -Ep). (39) 
The dimension of the stability matrix is (1 + 2N) 2 and it 

is sparse. For a given fixed point P A we compute numer- 
ically first few eigenvalues of the corresponding stability 
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TABLE I: Self-consistent electron densities on the molecule = P ss , (k = 1,2,3) for nonequilibrium fixed points computed 
by the Lindblad kinetic equation with different size of the buffer zones (N = 200, 400, 600) and by the NEGF method. 



u 


N = 200 


N = 400 


N = 600 


NEGF 


ni 


n 2 


n 3 


m 


n 2 




m 


n 2 


n 3 


m 


n 2 


ri3 


2.0 


0.34 


0.60 


0.64 


0.34 


0.59 


0.65 


0.34 


0.59 


0.65 


0.33 


0.58 


0.66 


2.05 


0.34 


0.54 


0.66 


0.33 


0.54 


0.67 


0.33 


0.54 


0.67 


0.33 


0.54 


0.67 


2.10 


0.34 


0.50 


0.66 


0.33 


0.50 


0.67 


0.33 


0.50 


0.67 


0.33 


0.50 


0.67 


2.15 


0.34 


0.47 


0.66 


0.33 


0.47 


0.66 


0.33 


0.47 


0.67 


0.33 


0.47 


0.67 


2.20 


0.34 


0.43 


0.66 


0.34 


0.44 


0.66 


0.33 


0.44 


0.66 


0.33 


0.44 


0.66 


2.25 


0.35 


0.40 


0.65 


0.34 


0.41 


0.65 


0.34 


0.41 


0.66 


0.33 


0.41 


0.66 




1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 
U 

FIG. 3: The dependence of Re(Ai) as a function of U for the 
unstable nonequilibrium fixed point. 



matrix with the largest real parts (Re(Ai) > Re(A2) > 
. . .). For this aim we used the ARPACK nonsymmeteric 
sparse eigenvalue solver. [l8j We find that the stability 

matrices resulting from pjo and P^l do not have eigen- 
values with positive real parts, therefore the correspond- 
ing fixed points are stable. Contrary, the stability matrix 

(2) 

resulting from P Q/3 always has at least one eigenvalue Ai 
with a positive real part, i.e., the corresponding fixed 
point is unstable. 

In Fig. [3] we plot Re(Ai) as a function of U for the un- 

— (2) 

stable fixed point P a g. The obtained curve is symmet- 
ric with respect to U = 2.143 where Re(Ai) reaches its 
maximum value. At U = 1.988 and U = 2.296, when the 
saddle-node bifurcation points are reached, eigenvalue Ai 
becomes zero. Thus, our stability analysis is consistent 
with the dynamical observations of Q that the "middle" 
fixed point can not be reached by time-dependent prop- 



agation of the NEGF Kadanoff-Baum equation. We do 
not observe purely imaginary eigenvalues of the stability 
matrix, therefore we rule out the existence of periodically 
oscillating fixed points in non-equilibrium Hartree-Fock 
approximation for the Anderson model. 

In Table IIIII we compare the calculated (N — 
200, 400, 600) fixed point densities of the molecule with 
the exact Hartree-Fock densities obtained by the NEGF 
method. The latter are the solutions of the following 
equation: [l9| 



1 f, T L (Lj)f L {u) + r R (w)f R (w) 

n = — I auj- 

it J (u> - 



Un - A(w)) 2 + (r(w)) 2 



(40) 



Here /l,_r(w) 



[1 



,(u>-nt,n)/T] 



is the Fermi-Dirac 



electron distribution in the left and right electrodes, and 
A = Al + An, T = + are the real and imagi- 
nary parts of the leads self-energy. For the tight binding 
electrodes, that we consider, the self-energy is given by 



Sl,h(w) 



XL 

2V? 



A i:fl (w) 




H 






-4V C 2 


ul,r + y 








Ak 2 - 





(41) 



where u>l,r = w — (-l,r- As seen from the table the 
larger the better our method reproduce the exact re- 
sult. For A^ = 400 our method quite well reproduce the 
exact result (the maximum deviation from the exact re- 
sults is under 2%), therefore our choice of the size of the 
buffer zone is justified and makes our Lindblad-type ki- 
netic equation numerically exact for electron transport 
calculations. 



J 



IV. CONCLUSIONS 



We have presented a theoretical method to perform a 
stability analysis of multiple nonequilibrium fixed points, 



which appear in self-consistent electronic transport calcu- 
lations. We employed the Lindblad kinetic equation and 
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underlying approximations of this kinetic equation were 
alleviated by the use of explicit buffer zone between the 
molecule and electronic electrodes. Nonequilibrium fixed 
point were obtained from the self-consistent solution of 
the kinetic equation in stationary limit. The asymptotic 
stability of these fixed points were studied by linearizing 
of the nonlinear kinetic equation. We obtained the non- 
Hermitian stability matrix from the linearized kinetic 
equation and analyzed its spectrum. If real parts of all 
eigenvalues of the stability matrix are negative, then the 
fixed point is asymptotically stable as t — > oo and it can 
be regarded as steady state. If at least one eigenvalue has 
positive real part, the solution becomes unstable and can 
not be the steady state. We obtained the explicit form of 
the stability matrix in Hartree-Fock and adiabatic time- 
dependent DFT approximations. The method was ap- 
plied to out of equilibrium Anderson model which yields 
three nonequilibrium fixed points under certain choice of 
parameters in nonequilibrium Hartree-Fock approxima- 
tion. We performed the stability analyse of these fixed 
points and demonstrated that one fixed point is asymp- 
totically unstable whereas the other fixed points corre- 
spond to physical steady states. 
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Appendix A: Lindblad kinetic equation for 
embedded molecule 



e iht ve -iht_ Now 



where Xl(t) = e x(i)e ' and vi(t) 
we introduce the density matrix for embedded system 
(molecule plus buffer) by tracing out the environment 
degrees of freedom 



p(t) = Tr EX (t) 



(A5) 



We assume that the total density matrix can be factor- 
ized and environmental degrees of freedom propagate in 
time as there were no interaction with the buffer (Born 
approximation) : 



XI (t) = Pl(t)p E - 



(A6) 



The density matrix for the environment is taken in the 
grand canonical ensemble form 



Pe ~ e "• keL e "• heR 

(A7) 

Setting /ix / piR and/or (3l ^ /3r in the environment, 
we drive the system out of equilibrium. 



After tracing, Eq. (|A4j) becomes 



Pi{t) = drTr E [[vj(t), [vi(t -r),pi(t- r)p E ] 
Jo 

Here, we have assumed that Tr E[vi(t), X/(0)] = 0- 

Using the explicit expression for the Hamiltonian h one 
can easily demonstrate that 



(A8) 



Let us begin with the Liouville equation for the total 
density matrix x(t) 

i*(*) = [W, *(*)]• (AI) 
We partition the Hamiltonian (J7J into two parts: 

h = H M +H e + H b + H sb (A2) 



v = H 



EB 



(A3) 



In the interaction representation the Liouville equation 
becomes 



iht — iht —is^t 

e afc CT e = e aka, 



e iht a ba e- iht = e~ lebt 



a ba + 0(l/N B ). (A9) 



where N B is the number of discrete single particle levels 
(i.e. the size) of the buffer zone. Therefore by choosing 
the large enough buffer zone, we may assume that buffer 
single-particle states evolve in time as free states. As a 
result, vi(t) takes the form 



vi(t) = ^2{vbk{t)al a a k<y + h.c. 

abk 



(AI0) 



1 



Mt) = 7 [vi(t), X im 



dr[[ W/ (i),b/(i-r), X/ (t-r)]], (A4) 



where Vbk(t) = Vbke^ 6k £i> '*. 

Now, the kinetic equation for the embedded system 
density matrix becomes 
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Pi{t) = - I dT Yl {^(*) u 6'fc(* - T ) C 1 - fkjal^avoiPiit - t)} - fk[a Va ,pi{t ~ T)]a\ a 

Jo akbb' 

: k {t)v Vk (t - t) f k a ba [al a ,pi(t - t)} - (1 - h)[a\, a ,pi{t ~ T)]a ba } 



r 



(All) 



Here f keL/R = Tr E (p E al a a ka ) = [1 + ^/»(«»-m/«3]-i 
Assuming that the environment relaxation time is very 
fast we can extend the integration range to +oo and 
p(t — r) ~ p(t) (Markov approximation). Finally, in the 
rotating wave approximation rapidly oscillating terms 
proportional to exp[i(e;, — e' b )t] for £t / e' b are neglected. 
Then, the kinetic equation (|Allj) becomes the standard 
Lindblad type master equation (ITo) . 

The obtained Lindblad master equation describes the 
time evolution of the open embedded system preserv- 
ing the probability and the positivity of the density ma- 
trix. Open boundary conditions are introduced via non- 
Hermitian dissipative part of Eq. (fT6]) , Tlp(t), which rep- 
resents the influence of environment on the system. The 
applied bias potential enters into Eq. (|16j) via fermionic 
occupation numbers /b (b e L, R) which depend on the 
chemical potential in the light and right electrodes. 



Appendix B: Stability matrix for nonequilibrium 
self-consistent DFT electron transport calculations 

Leaving apart the conceptual questions about the use 
of ground state DFT for self-consistent electronic trans- 
port calculations, we can say that the only difference 
in practical computations between nonequilibrium time- 
dependent Hartree-Fock discussion (Section II) and adia- 
batic time-dependent Kohn-Sham DFT is that the time- 
dependent Fock matrix (|24|) becomes 



F nm (t) =T nm + / dr0„(r)i£ s (r,t)0 m (r). (Bl) 



Here 



<•;.(.•./)= I dr'^±+v° c [p}(r,t), (B2) 



where v% c [p](r,t) is the exchange-correlation potential. 
Then the system of equations (|2T[) and (|29|) remains the 
same, but the steady state Fock matrix ([50]) becomes 

F nm = T nm + J dr<f> n {r)v% s [p](r)<f> m (r) (B3) 

and variation (I3T1) is 



SF nm (t) = J dr<j> n {r)6v a KS (r,t)4> m (r), (B4) 

where Sv^s ^ s gi ven by the standard expression of linear 
response time-dependent DFT 



dr 



,8p{v',t) 



J dv'f xc {av,a'v')5p al {v',t). (B5) 



The exchange-correlation response kernel / kc (ot,(tV) 
is given in the usual adiabatic approximation. [20l |21| 
i.e., the exchange-correlation contribution is taken to be 
simply the second derivative of the static ground state 
exchange-correlation energy E xc with respect to the fixed 
point spin density ~p a (r): 



f xc (ar,a'r')=S(v-r') — 



5E xc [p] 



5p a ,(r')5p a (v) 



(B6) 



Therefore, the stability analysis of nonequilibrium fixed 
points for self-consistent DFT electronic transport calcu- 
lations can be readily performed within standard adia- 
batic time-dependent density functional response theory. 

MM 
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